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ABSTRACT 

The occurrence of independent events at random in the 
plane, i.e. the formation of a planar point process, is 
discussed. Both homogeneous and nonhomogeneous processes 
are considered. A specific functional form for the parameter 
in a nonhomogeneous planar Poisson process is used to 
illustrate the development of test and parameter estimation 
techniques. The problem finds application in the description 
of biological phenomena as well as in search and detection 

problems . 
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I. INTRODUCTION 



Many problems arising naturally in a physical sense are 
often so complex that the identification and description of 
underlying mechanisms must use the tools of probability and 
statistics. Some of the reasons leading to the requirement 
of using these tools are: 

(i) the data base may be so large or complex as to 
preclude identification of any driving mechanism 
without recourse to statistical analysis; 

(ii) if identifiable, the mechanisms may be Inherently 
probabilistic; or 

(ill) if identifiable and deterministic, the governing 
law which the mechanisms obey may be unknown. 

This paper is concerned with the use of statistics in 
the identification and mathematical description of the spc ial 
distribution of events (occurrences). Included is the det c- 
tion and estimation of parameters which influence the 
description of this distribution. 

The area of concern here is a departure from those sta- 
tistical methods which have been developed to detect the 
effect of varying a controlled segment of the underlying 
mechanism. Among those methods would be the design of exper- 
iments, regression analysis, time series analysis, and 
analysis of variance. One goal of such analysis is to 
hopefully predict the advislbility of pursuing some course 
of action. 
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In the basic model of this paper, events are considered 
to occur vjith a Poisson distribution in the plane. This 
"is the natural model for the expression that 'points are 
distributed at random’," [Fisher, 1972, p. l4l]. The bi- 
variate Poisson process will be defined and then developed 
through the use of partial differential-difference equations, 
a widely repeated procedure in the univariate case but 
neglected in the bivariate case. 

Initially a homogeneous Poisson process will be assumed 
to control the underlying mechanisms. Then trends will be 
introduced by defining the Poisson parameter in such a way 
as to make it be spatially dependent. This will be the basis 
for the definition of the non-homogeneous Poisson process. 

Time inhomogeneity will not be considered. Thus, the data 
are assumed to be taken concurrently, i.e., the period of 
observation is short compared to any period of change of 
the parameters. 

Tests will be developed to distinguish between homogeneity 
and non-homogeneity and the method of maximum likelihood 
will be used to develop estimates of the parameter in the 
homogeneous case and parameters in the non-homogeneous case. 

In the latter case, conditional likelihood techniques will 
be utilized to develop tests and estimates. Throughout, 
testing and estimation procedures will be based on a single 
realization of the process which consists of the num.ber of 
events observed and their spatial locations. 
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The problem of concern finds application in the estima- 
tion of the density of trees in a forest; here one might be 
concerned with estimating the potential yield of lumber from 
a given forest area where inhomogeneities arise due to soil, 
weather patterns, topography and other physical reasons. 

Another application might be in naval search and detec- 
tion problems. For example, one might be searching for a 
merchant ship in distress whose location is not known exactly 
due to failure of the ship's communication equipment. Here 
the independence assumptions of the planar Poisson process 
may be valid, but not the assumption of homogeneity. In- 
homogeneities of location occur because of preferred sea 
lanes and physical characteristics of the ocean and 
atmosphere . 
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II. THE HOMOGENEOUS POISSON PROCESS IN THE PLANE (HPPP) 



A. GENERAL DE'^/ELOPMENT 

Consider a stochastic process of events occurring in 
the plane (i.e., a so-called point process) v/hlch is 
characterized by the assumptions 

I. There exists a finite positive constant X > 0. 

II. For any integer k ^ 1 and any set of non-overlapping 

regions R^ with areas A^,***,Aj^, (in the 

usual geometric sense) , the number of events occurring 
in any region i, denoted N(R^), has a Poisson dis- 
tribution with parameter XA^ which depends only on 
the area of the region, A^, and not its shape. Thus, 

( XA. ) ^exp(-XA. ) 

prob {N(R. ) = n. } = — — j — . (1) 

1 X i * 



III. 



Further, N(R^. ), i = l,2,*‘‘,k, are mutually indepe - 

dent in that N(R^) is not affected by the occurrer e 

of events in any other region or in any grouping of 

the regions, G, as long as R.Og = 0. Thus 

J 



prob{N(R^)=n^, 



l=l,--*,k}= 



k 

n 

1=1 



(XA^) 





( 2 ) 



Definition 1 : If a process obeys the above assumptions it is 

called a homo geneous planar Poisson process (HPPP). 

For reasons of arbitrary shape the above basic definitions 
will suffice. However, under certain geometrical assumptions, an 
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equivalent definition for the HPPP can be achieved in a man- 
ner similar to the development of the univariate Poisson 
process through the use of partial dif f erentlal-dlf feronce 
equations. This is useful for the development of statisti- 
cal properties and will be very Important in the development 
of the non-homogeneous process. Such a development: ar so 
provides another phenomenological approach to the ho.:, genecus 
Poisson process, one which might arise through the struc- 
turing of a model for instance. For Illustrative purposes 
the following development will be accomplished using rectan- 
gular regions. Note that the developmient is very dependent 
on the geometry Involved; hence developments with other 
geometries (e.g. circular regions) must proceed so;;.e-v'.as 
differently . 

The underlying assumptions in the differential equation 
development will be 

I*. There exists a finite positive constant \ > 0. 

II'. For any region R* with Incremental area 5.A, inde- 
pendent of the shape of the region except possibly 

as noted above 

(a) prob {no event in R*} = 1 - XAA + o(AA), 

(b) prob (one event in R*} = XAA + o(AA), 

(c) prob (more than one event in R*} = o(AA), 

where "g(AA) is o(AA)" means lim. = 0, or 

AA->0 

specifically in rectangular regions the limit as Ax 

^ ( A A ) 

or Ay or both go to zero of ? . - is zero. 
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Ill' . 



The occurrence of e\:nt-G in R* is independent 
of the occurrence of events in any region r"*" 
where R*0 r^ = 0. 

It will be shov.T! that I'. 11* and III' imply and are 
implied by I, II and III so that the two sets of assumptions 
are equivalent and hence the incremental assumptions give 
rise to a HPPP. Clearly :! cjk; I' are the same, as are III 
and III'. Also II implies II' since by ‘(1) 

2 

(a) prob {N(R«) = 0} = = 1 - XAA + ^ (AA)^ - ... 

2 ^ 

^ 1 - XAxAy + ^ Ax“Ay^ - . . . 

= 1 - XAA + o(AA), 

with the definition of o(AA) given above. Also 

(b) prob {N(R*) = 1} = XAAe = XAA(1 - XAA + ...) 

= XAA + o(AA) 

/,.,\1 ^-XAA 

and (c) prob {N(R«) > 2} = E = o(AA). 

1=2 

The problem remaining in order to demonstrate equivalence 
between the two sets of assumptions is to show that II' 
implies II. 

Consider a region R bounded by the co-ordinate axes and 
lines X = X* and y = Y*, with area X*Y*. Now extend the 
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sldeci to X = X*+Ax and y = Y*+Ay (see Figure 1). Consider 
the probability of n events occurring in the extended 
region, R' = R + where: 

(a) R has area X*Y*, 

(b) R^ has area X*Ay, 

(c) R^ has area Y*Ax, 

(d) R^ has area AxAy; 

(a)-(d) imply R' has area X*Y* + X*Ay + Y*Ax + AxAy. 

The assumptions I', II’, and III’ imply 

prob {no event in R^} = 1 - XX*Ay + o(X*Ay), 

prob {one event in R^} = XX* Ay + o(X*Ay), (3) 

prob {more than one event in R^} = o(X*Ay); 

prob {no event in R^} = 1 - XY*Ax + o(Y*Ax), 

prob {one event in R 2 } = XY*Ax + o(Y*Ax), (4) 

prob {more than one event in R^} = o(Y*Ax); 

and 

prob {no event" in R^} = 1 - XAxAy + o(AxAy), 

prob {one event in R^} = XAxAy + o(AxAy), (5) 

prob {more than one event in R^} = o(AxAy). 

Moreover, statements (3)> (^)> and (5) are independent. 

It is noted that the above equations may have two 
different interpretations. For instance in (3), prob{one 
event in R^} = X*Ay + o(X*Ay) is interpreted to mean one event 
in a two-dimensional process with parameter X and area X*Ay. 
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Figure 1. The incremental increase of a 
rectangular region. 
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However, another interpretation would be to consider the one- 
dimensional (marginal) process of events projected onto the 
y-axis, in which case the parameter is XX* and the incremental 
interval has length Ay. 

For rotational convenience, let P (X*,Y*) denote the 
probability that n events occur in a region with area X*Y*. 

The differential-difference equations are v;ritten noting that 
n events may occur in an extended region by having n events 
in the unextended region and no events in the extension, 
n-1 events in the unextended region and one event in the 
extension, etc. Hence 



P (X*+Ax,Y*) = P (X*,Y*) • P (Ax,Y*) 
n ’ n ’ o’ 

+ P^_^(X*,Y*) • P^(Ax,Y*) 

+ P^_2(X*,Y«) • P2(Ax,Y*) -r ... 

= P (X*,Y*)[l-XY*Ax] + P , (X*Y*)[XY*Ax] + o(Y*Ax). 
n n-i 

( 6 ) 



Similarly, 

P (X*,Y*+Ay) = P (X*,Y*)[l-XX*Ay] + P^_. ( X* , Y* ) [ XX*Ay ] + 
n n i* -L 

+ o(X*Ay), (7) 



and 
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P^(X*+Ax,Y*+Ay) = P^(X«,Y*)[l-XY^Ax][l-\X«Ay][l-XAxAy] 

*1 i i 

?^_^(X^'‘,Y^O [XY* Ax(l-XX*Ay ) (1-XAxAy) 

+ XX^-'Ay (l-XY*Ax) (1-XAxAy ) 

+ XAxAy (l-XX*Ay ) (1 -XY*Ax) ] 

+ P^_ 2 (X-^,Y*)[XX*Ay . XY*Ax(l-XAxAy) 

+ XX*Ay (l-XY*Ax)XAxAy 
+ (l-XX*Ay ) XY*AxXAxAy ] 

+ P^_3(X*,Y*)[X^X*Y*Ax^Ay^] (8) 

+ o(AxAy) + o(X*Ay) + o(Y*Ax). 

Interpreting the above equations, the third term on the 
right hand side of (8), for example, states that there can 
be n events in the extended region R* if there are n-2 
events in R and exactly one event in each of any tv;o of t e 
added regions. That is, there can be two events in the 
added regions R^ , R^ and R^ if one occurs In each of two 
regions and none occurs in the third region, i.e., one in 
R^, one in R 2 and none in R^, etc. Collecting all terms of 
order o(AxAy), o(X*Ay) and o(Y*Ax), (8) reduces to 
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P^(X*+Ax,Y*+Ay) = Pj^(X*,Y*)[l-XY*Ax-XX*Ay-XAxAy+X^X*Y^AxAy] 

+ Pj^_i(X*,Y*)[XY*Ax-2X^X*Y*AxAy+XX*Ay+XAxAy] 

+ P^_p(X*,Y*)[x^X*Y*AxAy]+o(AxAy)+o(X*Ay)+o(Y^Ax) . 

^ ( 8 ') 

= P„(X*,Y*)[l-XY*Ax]+P , (X*,Y*)[XY*Ax]+o(Y*Ax) 

n n-1 

+ P^(X*,Y*)[l-XX^Ay]+P^_^(X*,Y*)[XX*Ay]+o(X*Ay) 

- P^(X*,Y*)-XP^(X*,Y*)AxAy+XPj^_^(X*,Y*)AxAy 

+ X^X*Y*[P^(X*,Y*)-2P^_^(X*,Y*) + P^_2(X*,Y«)]AxAy-"o(.<- 

Noting from equations (6) and (7) that the fir-sr "hree 
terms on the right-hand side of the above equation arc- 
Pn(X*+Ax,Y*) and the next three are P^(X*,Y*+Ay) ar.:. 
rewriting (8’), the result is 

P^(X*+Ax,Y*+Ay)=P^(X*+Ax,Y*)+P^(X* ,Y*+Ay )-P^(X^ ,Y^ ) 

- XP^(X»,Y*)AxAy+XP^_^(X^,Y*)AkAj (8”) 

+ X^X*Y*[P^(X*,Y*)-2P^_^(X*,Y*)+P^_2(X*,Y*]AxAy 

+ o(AxAy). 

The definition of the second partial derivative with 
respect to two variables is 
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a c 3F(x.y) 

ay^ax 



}=iim _ 

Ay-vO^y Ax->0 



_ lini — } «Yy 

AX-.0 



Hence, transposing the first three tcrir.s of equation (8") 
to the right hand side, c:lvi.dir,r AxAy and then taking 
the double limit resultr 1; 



a F^(X*,Y*) ^ _^p )H XP 

axay ^ ' 



-2P T (X*,Y’-)-i-P 
n-1 ’ n- 



X'SY») + X^X*Y*[P^(X*,Y*) 
'X*,Y«-)] (9) 



The solution to (9) (a p^r... cifferential-dif ference 
equation) is easily shew ■ tc. 



P (X*,Y*) = K(XX*Y«)"e e^f'-,\X»^Y*) , n = 2,3, (10a) 

^ ■ n! ■ 

where K is an arbitrary constant. Special considerations 
are needed for n ='0,1 since for these cases some of the 
terms in (8") and (9) are not defined. Rewriting (8") and 
(9) while concurrently eliminating the proper terms leads to 



P (X*,Y*) = 
n * 



K(XX*Y*)’^exo(-XX*Y*) 

n! 



n = 0,1,2, . . . (10b) 



Since P^(X*,Y*) is a probability statement and for any 
given region the number of events in that region must be 
some non-negative integer, the constant K is seen to be unity. 
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Hence (10b) is equivalent to (1) which v;as to be proven. 

Thus the two sets of assumptions imply the same things, 
namely that the number of events in a region has a Poisson 
distribution with parameter proportional to the area of the 
region and independent of its shape and the number and 
position of events outside the region. Note that the 
formulation excludes multiple events, i.e., the occurrence 
of two or more events at any point or on any line in a single 
added region such as in Figure 1. 

Also a similar derivation will go through for circular 
regions using polar co-ordinates, but there are differences 
in the special properties of the Poisson process as defined 
through assumptions I, IT and III in differently shaped 
regions. These are discussed below. The differences in 
the special properties of the non-homogeneous planar Poisscr. 
processes as they vary with different geometries are an 
essential elem.ent of the analysis of points (events) in 
the plane . 

B. TESTING DATA FOR HOMOGENEOUS PLANAR POISSON PROCESS (,HPPP) 

Given the occurrence and spatial location of n events in 
a rectangular region of area X*Y*, consider the problem of 
determining whether or not these points occur as realisations 
of the KPPP. Miles [1970, p. 89] has stated a consequence of 
definition 1 as 

Corollary . Assume a rectangular region R^ with area A^. 

Given N(R^) = n and 0 < the n points are Independently 

and uniformly distributed in R^. 
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Proof: Let = X*Y* where is a rectangular region 

bounded by the coordinate axes, x = X* and y - Y*'. Label 
the n given points in any convenient manner, s .g. . on the 
magnitude of the y-component. Let (x,y)^^^ denote the i^^ 
labelled event. Consider an incremental region v/ith area 
dxdy which has the property: prob {exactly one event in the 

incremental region of area dxdy} = (dx,dy) = Adxdy + o(dxdy). 
Consider now n incremental recangles dx^dy^, 1 = l,***,n, 
placed in Ignoring probabilities of o(dx^dy^), assump- 

tions I, II and III imply that the Joint probability that 
t h 

the i event falls in the Incremental rectangle, dx^dy^, 
i = l,***,n, and exactly n events occur altogether in X*Y* 
is given by 

Xdx^dy^ . . . Xdx^dy^exp{-XX*Y* } . 

Restating in terms of the density function, 

f{(x,y) , . . . , (x,y) ,n;X) = X’^exp{-XX*Y*} , 

where f{...> is the Joint density of (x,y)^^^, i = l,...,n, 

and the probability that the number of events in X*Y* is n. 

The exponential term in the above expressions is an approx- 

n 

imation to exp{-(XX*Y* - I Xdx.dy. )}, i.e., represents 

i-1 ^ ^ 

the probability of no events within the region X*Y* but 
outside the incremental regions containing each event. 

By conditioning on the occurrence of n events in the 
region which are distributed Poisson with parameter XX*Y*, 
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. . (x,y)^^^ |n;X} = 



X%xp(-XX*Y*) 

(XX*Y»)^exp(-XX^Y*) 



n > 1 



n! 



n! 



( 11 ) 



(X*Y*)’^ 



which is the joint distribution for n bivariate uniform 
random variables ordered on one of the random variables as 
is shown in Appendix A. Note also the Independence of the 
conditioned density from the parameter X, i.e. the random 
variable N is a sufficient statistic for X. 

As a consequence of the above corollary, it is apparent 
that if the points of the HPPP, conditioned on the num*ber 
of events observed to occur, are in fact ordered v;ith 
respect to the increasing magnitude of the y-component, then 
no "information" is available about the ordering of the x- 
components, i.e., each of the n! orderings that can be 
induced on the x*s by the orderings on the y's has probability 
1/n!. This is readily apparent since in the bivariate uniform 
case the two components were independently selected. Hence, 
if (x,y)^j^^ 'is determined by i.e. the points are 

labelled by the ordered y-component, then 



prob{Xj^ = X^j^} = i j = 1,2,. . . ,n 



where is the j^ magnitude, and 




nj ^ lj***jJ IjJ'^lj***^} 







(12) 
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Hence if the x-components of the points ordered on the 
y-components exhibit any natural ordering then the x- and 
y-compcnents have not been independently selected and the 
observed process cannot be a HPPP. This v/ill be the basis 
for many of the tests for a HPPP against a non-homogeneous 
planar Poisson process to be discussed later. 

Lemma : If the bivariate process is Poisson and the regions 

are rectangular, then the projections of the events onto 
the coordinate axes may be shown to be univariate Poisson. 



Proof: Consider the occurrence of events in a rectangular 

region of area X*Y*. Then by III the occurrence of an event 
in an incremen*Cc,l strip is independent of all occurrences 
outside the strip. Hence the projections onto the coordinate 
axes give rise to independent counts along the axes. 



P (x.Y*) 

n' - 



( XY^x)‘^exp(-XY*x) 
n ! 



n = 0,1,... 
0 < X < X* 



and 

p (X* y) = 0-X*y)^exp(-XX*y) 
n ^ n ! 



n = 0,1,... 
0 < y < Y* 



which gives the univariate Poisson distributions v/ith 
parameters XY* and XX* respectively. 

Note here the inherent dependence on the shape of the 
assumed regions. In using rectangular regions equal lengths 
in the marginals reflect equal areas in the bivariate 
distribution. 
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If the regions were circular then vertical projections 
onto the axes would represent decreasing area as the dis- 
tance from the origin Increased. Since the occurrence of 
events Is assumed to be proportional to the area projected, 
an actual HPPP would Induce a non-homogeneous process on 
the marginals due to the distortion In the mapping. >’or 
clarification, refer to Figure 2. Hov;ever, If the re'_,:.ons 
are circular then radial projections could be utilised sc 
that the event occurring at In Figure 2 Is repre- 

sented In the x-marglnal by an event at X 2 » To def Inc- 

equal area projections In this case the transformation 
2 

X X = x’ Is made, In which case a unit Increase in x' 
defines the addition of a unit amount of area to the ■' '^Icr. . 
For example. If a unit area is generated by a circle- 
radius r = 1, then the area enclosed In the ring of 
1 _< r £ yr Is the unit area, as Is the area in the ring 
^ etc. In general, yrT £ r £ /n+1 defines in 

polar coordinates a ring with unit area. 

Returning to the assumption of rectangular regions, 
three characteristics of the HPPP are now available which 
can be used as the basis for testing a sample for belonging 
to the HPPP description of events In a rectangular region R: 

(A) Independence of' the x-orderlng from the ordering 
on the y-components . 

(B) Univariate HPP (homogeneous Poisson process) In 
the x-marginal and, conditionally on n events In R, a 
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-f 
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Figure 2. Vertical and radial projections of an 
event to form the marginal process. Shaded regions 
represent the deviations of projected areas arising 
from the rectangular projection of circular areas. 
Thus, the shaded regions Indicate the degree of 
non-homogeneity induced by the mapping. 
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uniform distribution of the distances to events. 

(C) Univariate HPP in the y-marginal and, conditionally 
on n events in R, a uniform distribution of the distances 
to events. 

Property (A) can be tested against general alternatives 
using a rank correlation procedure (or Spearman's correla- 
tion, see Pearson and Hartley [I 966 , Table ^^]). Properties 
(3) and (C) can be tested by standard univariate methods 
as in Cox and Lewis [1966]. 

Note that in the above discussion the interest lies in 
the nature of the process rather than in specifically 
describing the process. Thus the determination of the 
param.eter X of the Poisson process is not a current objec- 
tive and it can be considered to be a "nuisance" parameter. 
Hence the conditioning argument above and the resulting 
independence of the tests from the value of the parameter 
are justifiable. 

Now let be the probability of a Type I error gener 
ated in testing for randomness, ag be the corresponding 
probability in testing for HPPP in the x-marginal, and 
likewise for the y-marginal. Then the probability of not 
falsely rejecting the HPPP hypothesis due to the randomness 
test is 1 - a^, etc. Hence the combined probability of not 
falsely rejecting HPPP is 1 - prob {type I error} or 

1 - P(I) = (1 - a^)(l - ag)(l - a^). 
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Therefore 



P(I) = 1 - (1 - a^)(l - ag)(l - a^) 



(13) 



is the probability of falsely rejecting a HPPP hypothesis. 

If through physical considerations one of the tests seems 
more or less significant than the others, the analyst can 
choose the weightings to so reflect the physical properties. 
Otherwise the values (and thus the tests) can be vreighted 
equally. This need for the determination of weightings is 
the inherent disadvantage of a multi-level test. 

The individual tests proposed above will be briefly 
described. For the rank correlation test, consider each 
from (^>y)(j^) which is ordered on the y-component. Also 
consider the ordered realizations along the x axis, where 
x^ = ^ . Then 

6 E (i - (j),)^ 
i=l - 



r 



s 



1 - 




( 1^0 



where (j)^ is the position of in the x-ordered sequence, 
is the rank correlation statistic. 

The exact distribution for r can be approximated by 

o 

fitting a distribution to its moments as discussed by 

Kendall and Stuart [1951, p.^77]. The exact distribution 

of r is tabulated in Biometrika Tables For Statisticians 
s 

[ 1966 , Table 44, p.23] for observed values of n between 4 
and 10, and the introduction to these tables gives approx- 
imations for 10 < n < 20 and for n > 20. For 10 < n < 20, 
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can be treated as a product-moment correlation coefficient 

between normally distributed random variables. For n > 20, 

yn-1 is assumed to be unit normal. 

In testing the marginal distribution for HPP, two 

separate tests are proposed. First, the uniform conditional 

test is used to test against trends in the data. As stated 

in Cox and Lewis [1966, p. 153]> "IP the series has been 

observed for a fixed time t^ {length X*} and n events occur 

in (o ,t^) { (o ,X* ) } , then the uniform conditional test is 

based on the variables U/ . = T./t {= X/ . >./X*} (1=1 , . . . n) 

{i; 1 o {i; 

conditionally on N, being equal to n." The {brackets} are 

^o 

supplied to relate the material in Cox and Lewis [1966] to 

this specific problem, and N, = n means the number of 

^o 

occurrences observed is n. Note that in the conditioning — 
of the realizations the "nuisance" parameters XX* and XY* 
are eliminated. 

Secondly, a test based on the ordered inter-eveni; spac- 
ings is used to test Poisson against stationary event 
processes which may be non-Poisson. For this test, Durbin's 
modifications of the uniform conditional test is used [Cox 
& Lewis, 1966, p. 155]. Referring to Figure 3, Durbin's 
modification describes a transformation from the random 
variable X to the random variable T and then to the random 
variable S. 

Let = X* - If the X^^^, 1 = l,2,...,n, 

describe the "times to events" in a Poisson process, then 
the T^, i = l,2,...,n+l, are independent exponentially 
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Figure 3- The generation of the transformed variables 
from the original process 
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distributed random variables with parameter X. If the 
' s are then ordered and the S^'s are generated as shovm, 
then the S^*s are independent exponential random variables, 
v/here has the expectation l/( (n+2-1 ) X ) . Also the trans- 

I 

formation = (n+2-1 )S^ defines independent identically 

distributed exponential random variables with parameter 

. i t 

X, and therefore X. = Z S. , i = l,2 ,...,n defines the 

^ J=1 ^ 

times to events in a Poisson process with parameter X, 

tin, 

and U. = Z S. is the statistic upon which a new 

1 X J 

uniform conditional test is based. 

Both tests should be utilized as the uniform conditional 
test is more powerful when testing for trends while Durbin's 
modification is relatively more powerful in testing against 
stationary event process alternatives. However, these 
tests are not independent of each other and thus cannot 
be combined as in (13)* 

As an alternative to the above procedure, the region 

of concern may be partitioned into several sub-regions and 

the number of events in each subregion used as a basis for 
2 

X testing. This method is discussed by Kendall and Stuart 
[1951, pp. 57^-5] who mention the problem of choosing the 
"right" partition, adding "Whether a particular partition 
has statistical interest depends on the purpose of the 
analysis". Due to the underlying uniformity of the condi- 
tional distribution, this problem reduces to the selection 
of the number of regions which are then used to form equal 
area sub-regions. 
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Another alternative to the above testing procedure is 
the evaluation of the sample product -moment correlation 
coefficient under the bivariate uniform distribution. The 
procedure is discussed by Kov;alski [1972], but unfortunately 
his discussion does not address the b'’variate uniform 
distribution. Kowalski makes tv.’o points very strongly: 
"Firstly, the distribution of r (the sample product-moment 
correlation coefficient undei ncn-normial assumptions) may 
differ from its normal-theory form, and, secondly, v/e may be 
in a situation in which p is a poor measure of association." 
Hence, if the exact distribution for r under the bivariate 
uniform distribution were knovai, then an exact test for the 
HPPP (given the occurrence of n events) could be devised. 

Durbin [1970] has alsc’ uroposed distance methods for 
testing bivariate distributions. The process herein described 
is well-suited to the methods Durbin uses since he first 
transforms the observations so that they occur uniformly 
on the unit square. Hence the natural transformation 
x* = x/X^ and y' = y/Y* avoids the problem of possible lack 
of uniqueness which is the central objection to the use of 
distance methods. These methods allow the analyst to adopt 
Durbin’s bivariate analog of the Kolmogorov-Smirnov tests. 

The advantage of this method is the elimination of diffi- 
culties concerning multi-level tests and partitioning 
tests . 
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The tests described in this section are very general in 
nature, i.e. 

: the process is HPPP is tested against 

: the process is not HPPP. 

Hence the alternatives being tested against are multitudinous. 
If it is desired to test a realization as being from a HPPP 
against a specific form of departure from HPPP, better tests 
may be defined based on the nature of the specific alterna- 
tive. For instance, one such departure could be non-homoge- 
neity, i.e., where X is not considered to be constant but 
rather a function of location; this subject is considered 
in chapter III. Another departure might be in the nature 
of the process itself. For example, events may occur 
according to some fixed plan in which case the process is 
deterministic and thus non-Poisson. A process may develop 
in v;hich the occurrence of an event prohibits the occurre ce 
of another event for some interval about itself, In which 
case events are not independent of other events and are 
thus non-Poisson. 

It must be remembered, however, that tests against 
specific alternatives may ignore some features that a more 
general test would detect and thus each individual specific 
test applies only to the specific form of departure being 
considered . 

Moreover, in all reasonable stationary alternatives, 
it does not seem possible to derive the likelihood function 
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of the observations. One thus cannot derive exact tests. 

For tests against specific alternatives based on distance 
methods, see Holgate [1972]. Tests based on spectra are 
discussed by Bartlett [1964]. 

C. SIMULATING A HPPP 

Suppose one were concerned with searching for submarines 
which are assumed to be dispersed in such a manner that the 
locations at any moment are generated by a HPPP. If one 
search procedure is to be selected from, many proposed search 
procedures, then a possible manner of comparing the effec- 
tiveness of the proposed procedures is to utilize each pro- 
cedure against several simulated dispersions.. In such a 
simulation, the only "variable" which would be of interest 
would be the procedures, so all variables such as detection 
and classification parameters, facilities available, etc., 
would remain constant. Another problem which might be 
considered would be the effect of the change of such param- 
eters on the search procedure selected (i.e., a sensitivity 
analysis of the procedure to assumed operating 
characteristics and facilities). 

By the initial remarks of Section B above and the 
statement of equation (12), several methods of artificially 
generating realizations of a HPPP can be determined. These 
methods may then be utilized to simulate the HPPP. 

Assume that the parameter XX*Y* is given. To select 
the number N of events to be observed in the region with 
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area X*Y^ , generate a random number U distributed uniformly 
on .[0,1]. Set N =-• n if 



n-1 

E 

i 



dx^iY^)^exp(-XX*Y«) ^ (XX*Y*)^exp(-XX»Y*) 

1! < ^ 11 

^ (15) 



2 

The summations can be evaluated using either x or Gamma 
Integral Tables [Cox and Lewis, 1966, p.2^]: 



n-1 

Z 

1 



1! 



Prob {x2n > 2p} 



y 



00 n-1 -V 

/• V e , 

^ (n-1)! «''• 



Next, consider a random variable X distributed uniformly 
over (0,X*), denoted X U(0,X*), and another independent 
random variable Y U(0,Y*). As realizations of each 
random variable are generated, number them chronologically, 
i.e. in order of appearance. Generating n (as determined 
above) such realizations of each random variable yields 2n 
numbers; x^, . . . ,x^,y^, . . . ,y^. 

The final problem remaining is to select a scheme for 
mating the x- and y- realizations to form ordered pairs 
which will constitute the realization of the HPPP. A few 
such schemes are enumerated: 



1 . 



2 . 



The sequence forms a HPPP. 

If the y^ are ordered to form then the 

sequence forms a HPPP. 
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3. Similarly, forms a HPPP. 

Additionally, any random permutation of the in 2, 
the y^ in 3 or either random variable in 1 can be used 
to form a HPPP. Thus <(^n+l-l’^(i) forms a HPPP, etc. 

The goal of the simulation and the purpose of sim.ulating 
the process as a part of the overall analysis must nov; be 
considered. If during the simulation it is desired to 
generate independent realizations of the process, then each 
iteration must involve a selection of n, the drawing of 2n 
uniform variates and the mating of the variates through some 
scheme such as those outlined in steps 1-^J above. On the 
other hand, if it is desired to utilize variance reduction 
techniques, then for any drawing of 2n random variates 
several schemes could be used for the mating process. Here 
Independence is lost immediately and this loss must be 
balanced by some gain elsewhere in the analysis. 

D. ESTIMATION AND TESTING FOR THE PARAMETER FROM A HOMOGE EOUS 
PLANAR POISSON PROCESS (HPPP) 

If the hypothesis that the process is HPPP with some 
unknown value of the parameter X is accepted, one m.lght 
like to obtain a point estimate or confidence interval 
estimate for X, or to test that the process has som.e given 
parameter X^. Note that the parameter X, which was considered 
to be a nuisance parameter in the previous section where 
the structural aspects of the process per se were tested, 
now specifies the process completely. 
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Since, as was seen in Section B, it is possible to 
set up the joint probability density function of the 
observations in a HPPP, point estimation of X can be 
based on the method of maximum likelihood. Note, hov;ever, 
that each observation consists of a single "look" at (or 
realization of) the process rather than n observations of 
a single random variable. Since it is a stochastic process 
the observations are not independent and identically dis- 
tributed. Hence the usual justifications for maximum 
likelihood procedures are not valid; see Brown [1972] for 
extensions of maximum likelihood theory of estimation to 
realizations of a Poisson process. 

Using the results of Brown [1972], suppose ‘that n HPPP 
events are observed to occur in a rectangular region of 
area X*Y*. From (11), for n ^ 0, 

L = f { (x,y) , . . . , (x,y ) ^^^ ,n;X} = x’^e”^^ ^ (l6) 

In L = nlnX - XX*Y*, (0 <_ X < «'. ) 

If n 7 ^ 0, this function is at X = 0 and X = “ and since 

= ^ _ x*Y*. the slope of the function decreases 
dX X ’ 

monot onically from “ to -X*Y*. Thus In L has a unique 

H 1 T 

maximum at the point where = 0. Setting this 

derivative equal to zero yields a unique maximum likelihood 
point estimate for X as 

^ “ PY*"* (n ^ 1) . (17) 
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where X Is unbiased (since E(X) = = X) and has variance 

X/X*Y*. Note that as the observed area X*Y* becomes large, 
the variance of the estimate becomes small; thus, by 
Chebyshev's Inequality [Lamperti, 1966, p. 20] 

P{|X - X| > a} < , (a > 0) 

a 

and as X*Y* -> 



P{|X ~ x| ^a} 0 



for all positive a and hence X converges to X in probability. 

/s 

The latter statement is equivalent to the assertion that X 
is a consistent estim.ator for X. Also since the variance 
of X is X/X*Y*, X has an estimated variance X/X*Y* = n/(X*Y*)^, 
and an estimated standard error of /n7x^Y*. 

If n = 0, the above method is not applicable. In this 
case, it might be preferable to give a confidence interval 
estimate for X. Specifically, a one-sided test alternative 
is used to generate a test for the assumed value X^^^^ using 
as an acceptance region only n = 0. Intuitively, X^^^^ will 
be sm^all enough so that 1 (i.e., the expected 

number of observed events is less than 1). The hypothesis 
to be tested is : X = X > Defining 

a level of significance a from (l6) by 

-X X*Y« 

prob{N = 0]X = = 1 - a = e , 
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the hypothesis is accepted at the level a. Conversely, 
for any given value of a, >^^^^1 determined by 



^null 



X 5 fy* = In (1 - a) 



_ - In (1 - a) 

null 



where the determined is the largest value of X 

that the test will accept at the a level, given that n = 0. 

Returning to the case of n ^ 1 , in order to test that 
the parameter of the process has some given value Xq, assume 
that n events from a HPPP are observed in a region of area 
The hypothesis to be tested is X = Xq against 

the tv;o-sided alternative : X 7^ Xq although one-sided 
alternatives can also be considered. Since N is a random 
variable taking on all nonnegative integer values with som-'- 
positive probability for any Xq, there is always some 
possibility of an observed value of the random variable N 
(the observation being denoted n) falling outside any finice 
range of values. Thus a region (n~,n'*’) must be specified 
such that if N lies in the region the hypothesis Hq is 
accepted; otherwise the hypothesis is rejected. The level 
a of the test is the probability, given X = Xq, that N 
falls outside the region (n ,n'*'). 

Since the test has been defined to be two-sided, the 
level is split into upper and lower levels a"*" and a” 
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so that a 



= a + a . The procedure must consider values of 
A < Xq as well as values of X > Xq. To proceed. It Is 
necessary to define 



P+(n‘^;XQ) = P{N > n'^IX = X.} = a'^ 



(18) 



1 

» (Xr,X*Y*)^e ^ 

= Z . — 



0 =n 



T ! 



and 



P_(n ;Xq) = P{N < n"|X = Xq} = a (19) 

n- (X,X*Y*)-5exp(-V*Y*^ 

= E rn . 

j=o 

Thus, for a given a"*", an n"*" may be determined such that the 
statement (18) just holds. Also, for a given a", a n~ may 
be determined such that (19) just holds. 

The null hypothesis Is accepted at the a level If the 
observed value of N falls between the two prescribed limits 
(n"^ > n”), v;here prob{N (n ,n"*^)} = a. Note that as 
stated, the result Is Indeterminate since a, once given, 
leads to many values for a"^ and a” = a - a"*" which satisfy 
the given a. The manner of selecting and a must be 
stated. Arbitrarily It may bedeslrable to set a"*” = a = a/2. 
Asymptotically, this choice of a symmetric acceptance region 
Is reasonable since as n Increases, the distribution of N 
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is approaching the (symmetric) normal distribution. The 
choice of equal a"*" and a may not be reasonable, however, 
for small since the Poisson distribution is 

positively skewed. 

The statement prob{N i (n",n’’’)|X = = a is the result 

of the test of the hypothesis Hq*. A = Aq at a given, fixed 
level a. It is this result from which one must usually 
draw conclusions regarding specification of the process. 

If the information thus available, l.e. is rejected 
or accepted at the pre-determined a level, is deemed 
insufficient for the purposes of a decision maker (for 
example) then another possibility is that the post-analysis 
information may be extended by determining for each obser- 
vation the exact a, a^, at which the hypothesis would have 
been rejected. The decision maker is then left v;ith the 
problem of the determination of his own level of significance, 
possibly based on his intuitive grasp of the problem and 
its significance in a larger frame of reference. Once he 
has determined his preferred significance level, the hypoth- 
esis is rejected or accepted at the specified level by 
comparison with a^. Thus the decision maker has gained 
some influence over the analysis but has had to pay with 
some time to reflect on the problem at hand. Alternatively, 
he can use Informally as a "goodness of fit" of the 
hypothesis . 

Using (l8) and (19), the significance test is defined 
conventionally [Cox and Lewis, 1966, p. 30] to be; the 
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hypothesis X = v;ould be accepted at the level of signi- 
ficance a in a two-sided equi-tailed test if the observed 
number of events, n, is such that n, v/hen used alternatively 
in (18) and (19) (i.e., is assumed to be one or the other 
of the end-points of the acceptance region) , produces 
as a solution to 

P{n;XQ} = 2min{P^(n;XQ),P_(n;XQ)} = . (20) 

Note that each observed value of n generates a new for 
any assumed Xq , hence = ttg(n, Xq). For example, 

P(30;20) = .0^136, P(20;20) = .7628 and P(15;20) = .1332. 

It can be seen that 'the fixed level procedure is 
computationally simpler, since for a specified a and X^, the 
interval (n ,n"^) need only be computed once while in the 
latter procedure must be recomputed following each 
observation of N. 

The Inverse of the above approach v;hlch utilized the 

two-sided equi-tailed test of significance for a given 

value Xq leads to the determination of confidence interval 

estimates of X. Given that n events are observed, it is 

required to determine some limits on the range of X such that 
/ 

the true parameter value X* lies within the stated limits 
with a probability 1 - a. That is, it is required to 
establish a X”(N) and a X'*'(N) such that 

P{X"(N) < X* < X‘'’(N)|N = n} = 1 - a. (21) 
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Using P{N 1 n|X = x'*’} < 1 - a'*’ to define a x'*’ as the 
greatest X such that equality just holds and similarly 
using P{N £n|X=X}=a~to define a X~ establishes the 
limits such that (21) holds. For a proof of this, see 
Brownlee [1965, p. 121]. Note that for each realization 
of N, a new ordered pair (X jX"*") is defined so that the 
ordered pair is a function of a random variable and hence 
is itself a random interval. The procedure only states 
that for (1 - a) X 100^ of the observations the true 
parameter X* will lie v/ithin the limits selected. The 
limits for observed n from 0 to 50 are tabulated [Pearson 
and Hartley, 1966, Table 40]. 

For a normal approximation to the confidence interval, 
Cox and Lewis [1966, p. 31] define the upper a point of the 
unit normal distribution as c^, and give the relationship 



prob{-c^^ 

2 



< N - XX*Y* 




} = 



1 



a. 



the relationship being correct as XX*Y* -»• The confidence 

limits thus obtained are, to a second degree of approximation 
using a continuity correction and the estimate c(X) = yn/X*Y*, 



n 






+ 




For example, if 50 events are observed from a HPPP, the 
exact .05 confidence interval is 37.11 ^ XX*Y* ^ 65.92 
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[Pearson and Hartley, 1966, Table ^0] v;hereas the normal 

approximation gives 37*79 £ XX*Y* £ 66.07* 

2 

There also exist x approximations to the significance 
tests and confidence intervals [Cox and Leviis, 1966, p. 33 
Brownlee, 1965, p* 1733* 
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III. NON-HOMOGENEOUS PLANAR POISSON PROCESSES (IIKPPP) 



A. GENERAL DISCUSSION 

If the stochastic process described above is generalized 
to allov/ the probabilistic structure of the event process 
to be dependent on the location of the events, a non- 
homogeneous planar process is evidenced. In the simplest 
such case a non-homogeneous planar Poisson process (NHPPP) 
arises if, in the definition of the Poisson process given 
above, assumption I is modified to become 

I". There exists a positive finite function X(x,y) > 0. 
Also note that II is changed by the fact that the number of 
events in any region is not only a function of the area 
of the region, but also depends on the location of that 
region vflthin the universe under consideration. Thus X 
is now expressed as X " X(x,y), and assumption II becomes 
II". prob{N(R^) = n} 

exp{-A(A.)} 

nl 

where A(A.)= . / X(x,y) dxdy the symbol ./ Implying the 
1 A^ 

integral over an area and X(x,y) is assumed to be continuous 
over Rj. (with area A^) so that the integral is valid. 

Assumption III remains unmodified, i.e. events occur 
independently of any other event or collection of events. 

Under the additional assumption that X(x,y) is continuous 
within the region of consideration, the incremental 
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development of Chapter II may be extended to achieve a 
description of the NHPPP. Additionally the continuity 
assumption on X and the definition of the parameter in the 

f 

process as an Integral over X eliminates the difficulties 
of line discontinuities, although there may be cases v;here 
this is an important component of the problem. This problem 
is not considered here. 

Referring back to Figure 1 in Section II-A, consider 
specifically the Incremental strip defining region R^. If 
the strip is divided into n sub-regions of equal area by 
taking n equal Increments along the x direction each of 
length 6x, then, under the assumptions on the behavior or 
^(x,y), the process in the ith sub-region can be approximated 
by a HPPP with parameter X'= X(x,y)^ where (x,y)^ is an 
arbitrary point in the ith sub-region. Specifically (and 
arbitrarily) the lower left point is chosen for the 
succeeding discussion; thus the parameter for the first 
sub-region has parameter X = X(0,Y*). Continuing, the 
probability statements for occurrence of events become 

P 3 ^(x,y) = X(0,Y*)6xAy + o(6x,Ay) 0 < x <_ 6x, 

P^(x,y) = X(6x,Y*)6xAy + o(6x,Ay) 5x ^ x £ 25x, 



P]^(x,y) = X(j6x,Y*)6xAy + o(5x,Ay) j6x £ x £ (j+l)5x, 
where j = 0,1,..., n-1, Y* < y £ Y*+Ay and n6x = X*. 



Since the pi’obability of more than one event in 
is o(X*Ay) the probability statements above are additive 
and 

n-1 

prob {one event in R } = Z X( j 6x , Y* ) 6xAy + o(X*Ay). 

j = 0 

In the limit as 'n -»■ by the definition of an integral 

X5f 

prob (one event in R } = { / X(x,Y*)dx}Ay + o(X*Ay) . 

0 

By similar argument , 

Y* 

prob (one event in R„ } -- { / X(X*,y)dy} Ax + o(Y*Ax) 

0 

and 



prob (one event in R^} = X(X^' ,Y* )AxAy + o(AxAy). 

By comparison with equations (3), (^) and (5) the above 
statements lead to definitions for average parameters for 
each of the regions R^, R 2 and as 

X* 

A (X*,Y») = ^ / X(x,Y»)dx, 

X 0 

X,(X*,Y*) = ^ / X(X*,y)dy, ' (22) 

^ ^ 0 



and 



X^(X*,Y*) = X(X»,Y*). 



i|2 



Using these average parameters, the equations (3), (^) 
and (5) are generalized, resulting in the following 
statements : 

prob {one event in R^} = (X* , Y* )X*Ay + o(X*Ay), 

prob (one event in R 2 } = X 2 (X* ,Y* )Y*Ax + o(Y*Ax), (23) 

and 

prob {one event in R^} = X 2 (X*,Y*)AxAy + o(AxAy). 

Using the result (22) as defining the parameter in each 
of the incremental areas in Figure 1, equations (6), (7) 
and (8") become 



P (X*+Ax,Y*) = P (X*,Y*)[l-X„Y*Ax]+P ^ (X* , Y* ) X„Y*Ax 

n ’ n * 2 n-1 ’ 2 



+ o(Y*Ax) 



(2M) 




+ o(X*Ay) 



(25) 



and 



P^(X*+Ax,Y*+Ay) = P^(X*,Y*+Ay) + P^(X*+Ax,Y*) - P^(X*,Y«) 



- X„AxAy[P (X*,Y*) - P t(X*,Y*] (26) 

3 n n-1 



+ X^X2^*Y*AxAy[P^(X*,Y*)-2P^_^(X»,Y*) 



+ P^_2(X*,Y*)] 



+ o(Y*Ax) + o(X*Ay) + o(AxAy). 
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Rearranging terms, dividing by AxAy and taking the 
double limit as Ax -»■ 0 and Ay -»• 0 yields 

3^P (X*,Y») 

- A3[P^(X*,Y^) - P^_l(X^Y»0] (27) 

+ r3_X2X*Y*[P^(X»,Y«)-2P^_^(X*,Y*)+P^_2(X*,Y«)] 



which, together with the boundary condition that P^ is a 
probability statement, gives 



P 

n 



(X*,Y^) 



(A(X*,Y*))^ exp{-A(X^ ,Y*)} , n 
n! 



0 , 1 , 



( 28 ) 



where 



A(X^Y*) 

/ 



X* Y* 

/ / X(u,v) dudv. 

0 0 



(29) 



Thus, the number of events occurring in a region 
bounded by the coordinate axes, x = X* and y = Y* has a 
Poisson distribution with mean given by A(X*,Y*). The 
mean can be considered to reflect the cumulative effect of 
X(x,y) in the region of concern. 

If n events from a NHPPP are observed to occur in a 
rectangular region defined as usual with area X^Y*, and the 
events occur at (x,y)^^^, 1 = l,...,n, the labelling done 
on the magnitude of the y-component, then the joint density 
of the events and the probability that the number of events 
in X*Y* is n is given by 



( 30 ) 



n 

n X( (x,y) . . )exp{-A(X*,Y*) }. 
j=l 

Note that (30) is a direct generalization of (l6). 

Hence the NHPPP can be described in a fashion similar 
to the HPPP, but the expressions have acquired increased 
complexity due to the necessity for the inclusion of 
integrals to define the parameters. The degree of added 
complexity is dependent upon the choice of the specific 
functional form for X(x,y). The next section develops the 
expressions for one specific form. 

B. A SPECIAL CASE 

To consider the location dependent type of process, 
a particular form for X(x,y) is chosen as 

X(x,y) = exp {a+ Bx + yy + <Sxy}. (31'> 

Note that if Bx + yy + 6xy changes very little ovc-r the 
range of interest of x and y, then 

X(x,y) = (1 + Bx + yy + 5xy)exp{a). (32) 

Other relationships may be used; however, they may cause 
necessary and untidy restrictions on the values which the 
constants a, B, y and 6 may assume. In particular, X(x,y) 
must be greater than 0 and the bivariate exponential 
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polynomial (31) ensures this v;ith no restrictions on the 
range of the parameters. 

Additionally, algebraic manipulation of the form reveals 
that the curves of X(x,y) = c, c a constant, include some 
interesting properties. 

1. If 6 = 0, then In X(x,y) = c is a family of straight 
lines in the plane, intersecting the x-axis at an angle 
6 = tan ^(-B/y)* In this case a clock-wise rotation of 
the coordinate axes through an angle 0 would give an 
exponential function of y only. 

2. If 6 0, then In X(x,y) = c describes a system of 

contour lines which form a hyperbolic paraboloid v/ith a 
saddlepoint at (-y/ 6,-B/'S) as is Illustrated in Figure 
It may be helpful to interpret the Figure in terms of a 
section of forest which has been sampled. The line r 
describes a possible direction of steepest ascent (DSA) 
which passes through or near the region being sampled. 

This DSA may not be a topographic feature but rather a 
mathematical expression for a possible increase in density 
of trees along some line. Obviously, there may exist a 
strong correlation between this mathematical DSA and some 
topographic features. Note here that along the DSA maximal 
values for X(x,y) are found in the sense that departing the 
DSA at right angles leads to decreased values for X(x,y), 
i.e. decreases in the forest density. 
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Figure U. Contour lines for X(x,y) = exp{a+6x+Yy+<5xy } . 

Note asymptotes and region being described 
(hatched. ) Here 6/y = 2 , 6/6 = ^ and 

all coefficients are positive. 



( 



il7 



3. The exponential form can be extended v/ith little con- 
ceptual difficulty, but possibly greatly Increased mathe- 
matical difficulty, to describe a much v/lder range of 
possible circumstances. For Instance, It Is reasonable 

to assume that the DSA line will bend; hence terms such as 
2 2 

ex y and and higher order may be Included In the 

exponent . 

For the special form of (31), the cumulative or integrated 
intensity function A(X* ,Y* ) is given by equation (29) and 



where El (•) is the exponential Integral and Ei(*) = c + In ( • ) 



Jahnke and Emde [19^5, p. 2]. 

The likelihood function for the NHPPP may be develop ,*d 
in a manner similar to that used in the discussion of the 
HPPP. The discussion leading up to (l6) Is modified by the 
fact that the parameter is location dependent resulting in 



becomes 



A(X*,Y*) = 



^ A(X* ,Y*) . 5 
exp( a-By/'S)' 



^ E1{(y+ 5X*)(|+Y*) }-Ei{(Y+6x*)|} 

(33) 



- El{(3+6Y*)|} + El{^} 




n 



L = exp{-A(X*,Y*) } n X((x,y),.J, (n > 1) 

i = l ^ 



( 3 ^^) 



where (^>Y)(j_) ^ labelling of the coordinates of the n 

point events. Thus, 



In L = -A(X*Y*) + Z In X((x,y),.v) 

i=l 

follo^^s directly. For the special case of X(x,y) given by 
(31), equation (3^) becomes 

n n n 

In L = -A(X^'Y^')+3 I X. +Y 2 y,-+5 S x.y.+na (35) 

i=l ^ i=l ^ i=l ^ ^ 

where A(X^Y^) is given by (33). 

The above joint density, or likelihood, function pro- 
vides a functional form v;hich may be manipulated to accom- 
plish the two principal concerns of the analysis of point 
processes: hypothesis testing and parameter estimation. The 
obvious null hypothesis is : 3 = Y = iS = 0 , in which case 
the nonhomogeneous Poisson process is being tested for homo- 
geneity since a non-zero a yields a constant parameter 
X = exp(a) > 0. Should the above null hypothesis be rejected, 
then the analysis proceeds to develop estimates for the 
parameters 3, Y and 6. This phase of the analysis may 
proceed differently depending on hov'f many and which of the 
parameters were tested as being different from zero. The 
complete, and most complicated, situation develops when all 
parameters are determined to be non-zero. Testing of 
parameters is the topic of Chapter IV while Chapter V 
discusses the estimation of parameters determined to be 
non-zero as a result of the testing procedure. 
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IV. TESTING FOR NON-ZERO PARAMETERS 



A. PRELIMINARIES 

It Is desired to formulate a method for testing the 
data (l.e., the number of events and their locations) In 
order to determine which of the parameters In the model 
given by (35), specifically a, S, y and '6 , are non-zero. 

Note that three assumptions are Inherent at the outset: 
first, that the NHPPP model Is valid; second, that the 
testing for homogeneity In Section II-B led to the rejec- 
tion of the hypothesis of homogeneity; and third, that the 
physical phenomena can be modelled by the NHPPP given by 
(3^) with the parameter X(x,y) given by (31). 

Testing the Poisson hypothesis per se when the function 
X(x,y) Is not known Is a com.pound problem which will not be 
considered here. It is analagous to the compound problem In 
regression analysis of testing both for an unknown regression 
function and for independent equal variance errors. 

From the third assumption, the likelihood function for 
the data Is given by 

L = exp{-:;A.(X* ,Y*) expCBIx^. + yly^^ + 6Ex^y^} , (36) 

where, for clarity In the future development, C = exp{d}. 
Conditioning on the occurrence of n events, n ^ 1 and 
defining L{ (x ,y ) ^ ^^ , . . . , (x, y ) ^^^ | n ; X (x ,y ) } = L(n) leads to 
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L n! exp{S5:x^ + 

prob{N=n} ~ Y X 

iff exp{Bu + yv + 6uv} dudv)^ 

0 0 

(37) 



where L(n) is read "the likelihood function conditioned on 
the occurrence of n events." Note that conditioning on 
the number of events observed has resulted in an expression 
which is independent of the parameter ^(or a), i.e. for 
given 6, y, and 6, n is a sufficient statistic for a. This 
is convenient because a here is a "nuisance" parameter since 
the term.s of interest are those which would indicate non- 
homogeneity rather than the establishment of the overall 
rate of occurrence. Thus by using the conditional likeli- 
hood a may be eliminated and the testing can proceed for 
non-zero B, Y and 5. In other words the value of a should 
not influence the testing for non-homogeneity parameters. 

If n = 0, certainly no departure from homogeneity cou? i 
be evidenced and hence this case is covered by HPPP; see 
II-B above. Hence the case of interest is n ^ 1. 

Physically, the model (35) gives rise to a parameter 
surface X(x,y) which has the properties: 

(a) Bt^O; yj^O; InX forms a hyperbolic 

paraboloid superimposed on 
a tilted plane, i.e. some 
"warping" of the tilted 
plane is evidenced. 

(b) Bt^O; yj^O; 6 = 0: InX forms a plane, tilted 

with respect to the x-y plane. 

(c) 3=y = 0; 6?^0: InX forms a hyperbolic 

paraboloid . 
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(d) 6 = y='S = 0: 



In X forms a plane parallel 
to the x-y plane, i.e. a 
HPPP is evidenced. 



There are a number of possibilities for testing: 

(a) A test of 

Hq: 3=y=6=0 

against 

at least one of the parameters 3, Y, 6^0 

is a test for non-homogeneity which is more specific than 
those in Section II-B and is easily derived by likelihood 
ratio techniques. 

(b) The above test is not of great interest; generally the 
specific non-zero parameter is desired rather than just that 
at least one of the three is non-zero. This leads to the 
question of selecting the significant subset, a problem 
which is difficult and as yet is unresolved. 

(c) The simpler problem is to assume an ordering, i.e. that 
if 3 = Y = 0, the process is homogeneous (6 is then assumed 
to be 0) and if 3 or y is non-zero but 6=0, then higher 
order terms are assumed to be zero. However, if the test 
indicates non-zero 3 or y this may be due to an aliasing 
effect because of a non-zero 6. If further testing of 
5=0 against 6/0 reveals 6/0, then it may well be that 
the true situation is 3 = y = 0 but 6/0. The procedure 

to be followed will not discriminate this case. 
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The same aliasing effect occurs in testing of 6 = 0 
against 6^0 where g and y non-zero and it is desirable 
to perform this test without the effects of the non-zero 
B and y. These are thus nuisance parameters, as v;as the 
case with a in testing g and y. For the present model ( 35 ), 
one can eliminate these parameters because it is seen from 
the exponential form (36) that for any 6, (n, Zx^, Zy^) 
is a set of sufficient statistics for (aj g, y) . Thus 
6 = 0 is tested with some function of 2x^y^ given n, Zx^ 
and Zy^. This statistic has a distribution Independent 
of the parameters a, g, y. 

The reason for basing the conditional test on ^x^y^ 
is that this is (conditionally) a sufficient statistic 
for 6 . 

B. SPECIFIC TESTS 

Assuming that some ordering exists on the parameters as 
discussed in possibility (c) above, tests are performed 
using the sufficient statistics (n, Zx^, Zy^^, Zx^y^) to 
determine if any non-homogeneity is evidenced by the data 
(l.e., through the statistics). This testing is more 
specific in nature than the testing encountered in Section 
II-B above due to the selection of a particular model. 

The set of sufficient statistics arises from this choice of 
a specific model to use as an alternative to homogeneity. 
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The testing vrlll assume the follov/ing sequence: 

(1) Condition on n and set 6 = 0. Test B = y = 0 

against B 0 or y 5^ 0* Note that it v/ould not be 

informative to test either B or y as a separate entity since 
in the formulation of the model B and y are unique only up 
to an angle of rotation. That is, testing of B and y 
jointly amounts to the detection of any tilt in In X(x,y) 
with respect to the x-y plane, regardless of the direction 
of the tilt. Failure to reject leads to the assumption 

of homogeneity due to the assumed ordering. 

(li) Rejection of leads to testing of 

Ho(ii): 6 = 0, < B < “ and -oo < y < oo 

against 

^l(il)" -co< g <00 and -oo < y < oo. 

The test thus specifies y and B as nuisance paramieters. 

In this test it is necessary to first condition on n. Ex 
and Ey^ to eliminate the nuisance parameters. 

In (i), conditioning on n and setting 5=0 leads to 

(By)'^ n! expCBEx, + yEy, } 

L(n) = i . 

(exp{BX} -1)^ (expCyYl - 1)" 

From this it is seen that the statistics (Ex^,Ey^) are 
(conditionally) jointly sufficient statistics for B and y. 
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Under 

Ex^/n + N(X*/2,X*^/12n) and Ey,/n -► M(Y*/2,Y*^/a2n) 

and the statistics are independent (see Section II-B) . Hence 
the expression 

Zx /n - X*/2 2 Zy./n - Y*/2 2 

— ) + (-i-— — ) 

x*//i^ y*/ 

2 

is asymptotically X 2 * Rejection or acceptance of 

is based on the adherence of the calculated value of this 
2 

sum to the x distribution , i.e. Hq is accepted if this 
sum has sufficiently small values. Acceptance of Rq^)j 
as stated earlier, leads to assumption of HPPP; refer to 
Chapter II. 

Following the rejection of It is necessary to 

proceed with testing of RgCii)* seen from, an 

examination of (37), the complexity of the exact distribu- 
tion following another conditioning argument (i.e. on 
n, Ix^ and Zy^) is prohibitive. However, for large sample 
sizes the conditional distribution can be approximated from 
the fact that Zx^/n, Zy^^/n and Zx^y^^/n, conditioned on n, 
are Jointly normally distributed for large n. Thus the 
asymptotic distribution of Zx^y^/n, given n, Ix^/n and 
Zy^/n, can be found from normal theory multiple regression 
results . 
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Under the assumption that 3 = y ~ ^ Oj the trivariate 



normal distribution v;hich arises is characterized by a vector 
and a matrix. The vector (]^) of expected values and the 
variance-covariance matrix (^) are given by 




and 



E = 



^X^/12n 



,X^y/2i4n 



0 

Y^/12n 

XY^/24n 



X^Y/24n 

XY^/24n 

7X^Y^/l44n 



from which ~ P 13 “ P 23 “ 0 . 65 ^ 65 . 

In the model given above, 

^O(ii)' ^ -03 < 3 < 00 ; -CO < Y < 00 

is to be tested against 

^l(ii)* *57^0; _co<3<co; _oo<Y<«». 



Since 2 x^y^ is a sufficient statistic for 6 when n, Ex^ 
and Ey^ are given, the test can be based on Ex^y^. Its 
asymptotic (conditional) normal distribution has mean 



I 



V 



xy 



and standard deviation a 



xy 



given by 
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t 



y 



xy 




a 

^ (Ex./n - X/2) + 



+ 




(Zy^/n - Y/2) 



and 



I 2 

a =0 (1 -Pt_ 

xy xy ^13 



2 . 
P23 ) 



E X . y . /n - y ' 

Thus under ^ is distributed as a unit 

^ ^xy 

normal variate and is accepted if this statistic has 

sufficiently small values. Failure to reject would 

imply that the In l(x,y). plane is tilted with respect to 
the x-y plane, but no "warping" is evidenced. 

The above development relies heavily on asymptotic 
assumptions. Small sample problems will be much more 
difficult to analyze. Any point in the above procedure wJ ch 
lead to rejection of any hypothesis would require the anal 'sls 
to proceed with the estimation of the non-zero parameters 
This is the subject of the next chapter. 
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V. ESTIMATION OP PARAMETERS 



It is desired to formulate a method for estimating the 
parameters a, Bj y and 5 of the non-homogeneous planar 
Poisson model given in IV-A where it has been established 
that a non-homogeneous process is evidenced by the data. 

Taking the logarithm of the conditional likelihood 
function (37) results in 

In L(n) = In n! + BZx^ + yZy^ + 6Ex^yj^ + n In A, (38) 

where A =A(X*,Y*)/^. Point estimation of (a, 6, y, 6) 
by the method of maximum likelihood uses the conditional 
likelihood function (38) to develop the estimates. See 
Section II-D for comments regarding use of maximum likelihood 
in this application. The solution to the set of simultaneous 
equations 

Y* X* 

Zx. - T ^ / u exp(Bu + yv + 6uv} dudv = 0 

^ '^0 0 

Y* X* 

Zy . - T ^ / V exp{Bu + yv + 6uv} dudv = 0 (39) 

^ '^0 0 

Y* X* 

Zx.y. - r / uv exp(Bu + yv + 6uv} dudv = 0, 

1 1 A 0 0 ■* 

/\ 

if obtainable, provides the point estimates B, y and 5. 

Note that this approach neglects the homogeneous term 
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during the estimation of the parameters giving rise to 
non-homogeneity. The neglected parameter may be estimated 
last . 



In order for the solution (g, y> *5) to equations (39) 
to describe a relative maximum to In L|n, it is necessary 
and sufficient that the matrix of second partial derivatives 
(£) be negative definite, see Frisch [1966, p. 120]. In 
examining this matrix in the case of (38), it is helpful 
to define S(u,v) = exp {gu + yv + 6uv}. Then the function 



has the properties: 

(a) s(u,v) ^ 0 

y* x«- 

(b) / / s(u,v) dudv = 1 . 

0 0 

(c) s(u,v) is continuous on [0 u £ X*, 0 ^ v ^ Y*]. 

Hence s(u,v) is a probability density function [Gnedenko, 
1962, p. 171]. 

Hence the matrix ^ can be shown to have diagonal 
elements such as 



Y* 




a 



11 



= - n[ / 

0 0 



0 0 



n Var U. 
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Continuing, the result is (where V/ is defined to be the 
W = UV) 



E 



( 



-n Cov (U,V) 
-n Cov (U,W) 



-n Var U 



-n Cov (U,V) -n Cov (U,V/) 

-n Var V -n Cov (V,V/) 

-n Cov (V,V/) -n Var W 



and ^ is revealed to be a covariance matrix. Note that the 
condition for a relative maximum, i.e. ^ negative definite, 
is independent of the realizations. 



matrix for a tri-variate distribution. But ^ is positive 

semi-definite [Gnedenko, 1962, p. 212], hence is negative 

semi-definite. That each of the principal minors has 

non-zero determinants remains to be shown. 

By the expressions given in Gnedenko [1966, p. 212], 

the covariance matrix ^ can be seen to be a Kankel matrix 

[Gantmacher, 1 , 1959, P- 338]. Hence if the rows of £ ar' 

linearly independent, then the determinant of £ > 0. But 

also Var U > 0 since U is a random variable and Var U Var V - 
2 

Cov (U,V) > 0 since the case of line discontinuities has 
been excluded (i.e., U cannot be a linear function of V). 

By the same reasoning, W is linearly independent of U and 
V. Hence all principal minors are greater than zero, hence 
£ is positive definite, hence £ is negative definite. 

Thus (3, Y, '5) provides ^ least a relative maximum to 
In L I n. 



Now E = -n E where E is the usual variance-covariance 
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If it were possible to determine that (6, y, 6) provides 



a global maximum to In L|n in the region of Interest, then 
conclusions as to uniqueness of the estimator could be 
drawn. Unfortunately, global extrema are difficult to 
establish. Since the method of estimation used was maximum 
likelihood, the estimates are consistent. Questions of 
biasedness are unresolved. 

In order to solve the system of equations (39), it is 
necessary to determine initial values for the parameters as 
a starting point for an iterative procedure. The partial 
differentiation of InL (35) with respect to the parameters 
and setting these partlals equal to zero results, after 
some algebraic manipulation, in 



n -A(X,Y) 



0 




6 





ea ^^BX(^(y+5X)Y 

X I- r/XAIY 



6 Y+<SX 




Zx y + ( ) A t 
^ ^ 6 



6(6+6Y) (y+SX)"- 6 






- - 1] - 6X[e^^ - 1] - - 1]} =0 



If it is assumed that the sum gX + yY + 6XY is small 



(near zero) as well as the individual terms in the summation 
being small, then the exponentials can be approximated by 



exp{x} = 1 + X, X near zero. Using the first equation 
in system (40) to give the Value for A(^*' jY* ) , i.e. 

= n, and the linear approximation in the remaining 
terms gives the abbreviated system: 



SXi + ^ n =0 

Zyi + |- n =0 (4l) 



y . + n 

11 



(1+g) 

5 (e+<5YKy+6X) 



+ 66X^Y + y<SXY^ 
-6^XY + eyXY] = 0 



The solution to (4l) provides the initial estimates for the 
parameters. These estimates can then be used in (39) or 
(40) to search for sequentially closer and closer 
approximations in a mathematical programming approach. 

Following the determination of the estimates S, y 
and 6, 5 can be determined from the solution to the first 
equation in the set (40). 

The determination of confidence intervals and levels 
of significance is not considered. 
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VI. CONCLUSIONS 



The procedures in Chapters III-B, IV and V are dependent 
on the particular choice of parameter form; hovfever, vrith 
different forms the concept of a non-homogeneous planar 
Poisson process may be used to describe a v;ide variety of 
"randomly” occurring phenomena. The choice of parameters 
which may be used is limited only by assumption I, l.e. 
positivity. One advantage of the method discussed herein 
over previously proposed schemes is the fact that the • 
specific form used admits the possibility of a ridge or 
line of maximum density to be mathematically specified 
and estimated. 

Also there is an attempt to describe the underlying 
process that caused the points to appear where they did, 
as opposed to using, for instance, the arc within which 
the most events were observed as the point estimate for 
the direction of maxim.um increase. 

Further efforts in this area Include a generalization 
into four dimensions (x,y,z,t) in order that zoological 
as well as botanical densities may be studied. Of especial 
interest is the estimation of densities of aquatic life 
and how the observed density fluctuates with season and 
with changes in environment. The latter problem seem.s of 
prime importance in evaluating the effects of anti-pollution 
programs on the fluid systems in which plants and animals 
exist . 
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Another problem which is closely related to the above 
Is that of Imperfect sampling and how the estimates are 
biased by sampling techniques. 

Chapters III, IV and V may be redefined In terms of 
data -gathered within a circle about some fixed point, 
especially with consideration of the relative efficiency 
of this data form referred to by Matern [i960]. 
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APPENDIX A: THE BIVARIATE UNIFORM DISTRIBUTION 



2 

Given a region R in E of area A and the fact that the 
probability of occurrence of an event in any sub-region R^ 
of area A^ within R is simply A^/A, a bivariate uniform 
distribution is described. For definiteness assume the 
region R is rectangular, so A = X*Y*. Nov^ 

Prob (X £ x,Y < y) = ■ = Prob (X < x) Prob (Y < y) 

X*Y* 

for 0 ^ X £ X* , 0 £ y £ Y* , in which case it is apparent 
that the coordinate axes define Independently chosen 
univariate random variables. 

Also, the density function is immediately 

f(x,y) = 1/X*Y* 0 < X < X*, 0 < y < Y*. 

Prom the density function the joint density for n independent 
bivariate uniform random variables is 

f((x,y)^,. . . ,(x,y)^,n) = 1/(X*Y*)’^ 

t h 

where (x,y)^ denotes the i ' pair of random variables 
selected. Now n pairs of random variables, or more simply 
n points in the plane, can only be ordered (without 
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replacement) in n! ways. Independent of the ordering process 
chosen . Hence, the joint density function for n ordered 
bivariate uniform random variables is 

f((x,y)(^),...,(x,y)(^^,n) = n!/(X*Y*)'' 

t h 

where ^x,y)^^ is the i^ point selected in the ordering 
scheme utilized. 

As a specific example, consider the n points to be 
labelled with respect to increasing magnitude of the y- 
component. Then 

" ^(k) ^ and (x,y)^j^^ = 

If the x-components are also ordered, then the set of 

p 

points P = ( (x^ ,y ^ j ^ i,j = l,...,n) defines n“ points, 
of which n are known to-be "occupied," that is, to descri e 
an event. For there exists some j such that Yq- ) 

gives the y-coordinate value for the event which gave ri£.e 
to x^^^. Similarly, for x^2) there are 

now n -1 j's remaining, one of which must correspond to the 
event giving rise to x^2)* Continuing to x^^^, there can 
only be one J left to be associated with the last x-value. 
Thus there are n! combinations of (x,y)^^^, i = l,...,n 
each having density of l/CX^Y*)*^ and so the ordered bivariate 
uniform density is established as that stated above. 
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